here::i_am("beo_replication_code.R")
library(here)

##Reading in all necessary data

beo <- read.csv("./Data/beo_full.csv")

##Loading necessary libraries
library(lme4)
library(stargazer)
library(ggplot2)


##Table 1: Effect of Total Number of BEOs on Black Turnout

table1 <- lmer(black_vote_pct ~ total + total2 + pop_thousands + inc_thousands + 
                  college_ed + black_pct + urbanpct +TS1970 + pres_elec_year + (1 | name) , data=beo)

stargazer(table1, dep.var.labels = "Black turnout percent", covariate.labels = c("Number of BEOs", "(Number of BEOs)$^2$", "Population (thousands)",
                                                                                 "Median income (thousands)", "Percent with college degree", "Percent Black",
                                                                                 "Percent urban", "Year", "Presidential election"))
##Figure 2: Predicted Black turnout by number of BEOs
#Note: Figure 2 relies on bootstrapped confidence intervals that are significantly computationally intensive. Code for bootstrapping can be found in beo_bootstrap.R
numCIs <- read.csv("./Data/number_cis.csv")

lwr_loess <- loess(lwr ~ val, data=numCIs)
upr_loess <- loess(upr ~ val, data=numCIs)
smoothed_lwr <- predict(lwr_loess)
smoothed_upr <- predict(upr_loess)

pdf("./Figures/numb_beos.pdf")
ggplot(numCIs, aes(x=val))  + ylim(25,80) + xlim(0,45) + labs(x="Number of BEOs", y="Predicted Black Turnout") +
  geom_ribbon(aes(ymin = smoothed_lwr, ymax = smoothed_upr), colour="grey10", fill="grey1", alpha=0.1) +
  geom_line(aes(y=pred), color="black", size =1 ) +
  geom_rug(data=beo, aes(x=total), alpha=0.2, inherit.aes=FALSE) +
  theme_bw()
dev.off()

#Table 2: Effects by BEO Type on Black Turnout

beo$city_dummy[beo$city==0] <- 0
beo$city_dummy[beo$city>0] <- 1

beo$county_dummy[beo$county==0] <- 0
beo$county_dummy[beo$county>0] <- 1

beo$state_dummy[beo$state==0] <- 0
beo$state_dummy[beo$state>0] <- 1

beo$judicial_dummy[beo$judicial==0] <- 0
beo$judicial_dummy[beo$judicial>0] <- 1

beo$education_dummy[beo$education==0] <- 0
beo$education_dummy[beo$education>0] <- 1

table2 <- lmer(black_vote_pct ~ city_dummy + county_dummy + state_dummy + fed + judicial_dummy + education_dummy + pop_thousands + inc_thousands + 
                  college_ed + black_pct + urbanpct +TS1970 + pres_elec_year + (1|name), data=beo) 

stargazer(table2, covariate.labels=c("City BEO","County BEO","State BEO", "Federal BEO", "Judicial BEO", "Education BEO", 
                                      "Population (thousands)", "Median income (thousands)",
                                      "Percent with college degree", "Percent Black",  "Percent Urban", "Year", "Presidential election"))
beo$nob2 <- beo$imputed_nob^2

model21 <- lmer(black_vote_pct ~ city_dummy + county_dummy + state_dummy + fed + judicial_dummy + education_dummy + pop_thousands + inc_thousands + 
                  college_ed + black_pct + urbanpct +TS1970 + pres_elec_year + (1|name), data=beo) 

#Figure 3: Predicted Black turnout by proportion of BEOs
#Note: Figure 2 relies on bootstrapped confidence intervals that are significantly computationally intensive. Code for bootstrapping can be found in beo_bootstrap.R
simCIs <- read.csv("./Data/proportion_cis.csv")

lwr_loess_prob <- loess(lwr ~ val, data=simCIs)
upr_loess_prob <- loess(upr ~ val, data=simCIs)
smoothed_lwr_prob <- predict(lwr_loess_prob)
smoothed_upr_prob <- predict(upr_loess_prob)

pdf("./Figures/prop_beos_revised.pdf")
ggplot(simCIs, aes(x=val))  + ylim(25,80) + xlim(0,0.437) + labs(x="Proportion of BEOs", y="Predicted Black Turnout") +
  geom_ribbon(aes(ymin = smoothed_lwr_prob, ymax = smoothed_upr_prob), colour="grey10", fill="grey1", alpha=0.1) +
  geom_line(aes(y=pred), color="black", size =1 ) +
  geom_rug(data=beo, aes(x=total_proportion), alpha=0.2, inherit.aes=FALSE) +
  theme_bw()
dev.off()

#Table 3: Effect of Number of BEOs not on ballot on Black turnout

beo$city2 <- beo$city^2

table3 <- lmer(black_vote_pct ~ city + city2 + pop_thousands + inc_thousands + college_ed + 
                  black_pct + urbanpct + TS1970 + pres_elec_year + (1|name), data=beo)

stargazer(table3, covariate.labels=c("Number of off-ballot BEOs","(Number of off-ballot BEOs)^2","Population (thousands)", "Median income (thousands)",
                                      "Percent with college degree", "Percent Black",  "Percent Urban", "Year", "Presidential election"))

##End of main text replication code
